package rcs34.libs.Sound;

/**
 *
 * @author RCS-34
 *
 * This is the Cooley-Tukey Fast Fourier Transform algorithm implementation
 * class. It computes FFT same as Inverse FFT over given samples. It also provides a
 * simple interface for the frequencies denomination.
 * Class is designed to be quick enough to use for a real-time signal processing.
 *
 */
public class FFTCooleyTukey {

  private static final double PI2 = Math.PI * 2;
  private final int sampleSize;
  private final int sampleRate;
  private final double[] freqs;
  private final double[] sins;
  private final double[] coss;

  /**
   * Constructor of the FFT processor.
   *
   * @param _sampleSize How many samples it will work with. Note, that the
   * sampleSize in FFT co-determines the nominal frequencies spectrum.
   * @param _sampleRate Sample rate of the input signal. The processor needs to
   * know that to be able to provide the frequency-based functionality
   * @throws Exception The exception is threw in case that the _sampleSize is
   * not a power of 2
   */
  public FFTCooleyTukey(int _sampleSize, int _sampleRate) throws Exception {
    sampleSize = _sampleSize;
    sampleRate = _sampleRate;
    double log = Math.log(sampleSize) / Math.log(2);
    if (sampleSize != Math.pow(2, log)) {
      throw new Exception(_sampleSize + " is not power of 2");
    }
    //
    int cacheSize = sampleSize / 2;
    sins = new double[cacheSize];
    coss = new double[cacheSize];
    freqs = new double[cacheSize];
    for (int i = 0; i < cacheSize; i++) {
      double iPart = i / (double) sampleSize;
      freqs[i] = sampleRate * iPart;
      double angle = -PI2 * iPart;
      sins[i] = Math.sin(angle);
      coss[i] = Math.cos(angle);
    }
  }

  /**
   * Method that can reverse-compute the position of the given-frequency-related 
   * bin in the output array
   *
   * @param _freq The given frequency
   * @return returns the relative position of the given frequency in the output
   * array
   */
  public double getIndexOf(double _freq) {
    return (_freq * sampleSize) / sampleRate;
  }

  /**
   * Returns the nominal frequency on the given index
   *
   * @param _index Index in the output array
   * @return Frequency in Hz
   */
  public double getFreq(int _index) {
    return freqs[_index];
  }

  /**
   * Gives the sample size which this processor was created with
   *
   * @return Sample size
   */
  public int getSampleSize() {
    return sampleSize;
  }

  /**
   * Compute the FFT over the sample. It takes the input values from the _sample
   * array, computes the FFT and saves the result to the _result array. The both
   * arrays size have to be of size [>2][sampleSize].[0] holds the real parts,
   * [1] holds the imaginary parts.
   * 
   * @param _input Input sample (RO). In case, that the input data are not
   * complex (audio samples), the imaginary parts should be set to zero.
   * @param _output Result (RW). Result of the FFT is stored there.
   * @throws Exception If _sample of _result does not the size of sampleSize.
   */
  public void fft(double[][] _input, double[][] _output) throws Exception {
    if (checkParams(_input, _output)) {
      ditfft2(_input, 0, _output, 0, sampleSize, 1, 1);
    }
  }

  /**
   * Compute the Inverse FFT over the sample. It takes the input values from the
   * _sample array, computes the FFT and saves the result to the _result array.
   * The both arrays size have to be of size [>2][sampleSize].[0] holds the real
   * parts, [1] holds the imaginary parts.
   *
   * @param _input Input sample (RO).
   * @param _output Result (RW). Result of the Inverse FFT is stored there.
   * @throws Exception If _sample of _result does not the size of sampleSize.
   */
  public void ifft(double[][] _input, double[][] _output) throws Exception {
    if (checkParams(_input, _output)) {
      ditfft2(_input, 0, _output, 0, sampleSize, 1, -1);
    }
  }

  private boolean checkParams(double[][] _input, double[][] _output) throws Exception {
    if (_input.length < 0
            || _input[0].length != sampleSize
            || _output.length < 2
            || _output[0].length != sampleSize
            || _output[1].length != sampleSize) {
      throw new Exception(
              "Sample size and/or result size does not match the value1 "
              + sampleSize + " for which the FFT was constructed");
    }
    for (int i = 0; i < sampleSize; i++) {
      _output[0][i] = 0;
      _output[1][i] = 0;
    }
    return true;
  }

  private void ditfft2(double[][] _input, int _inOffset, double[][] _output,
          int _outOffset, int _size, int _stride,
          int _inversion) {
    if (_size == 1) {
      if (_inversion < 0) {
        _output[0][_outOffset] = _input[0][_inOffset] / sampleSize;
        _output[1][_outOffset] = _input[1][_inOffset] / sampleSize;
      } else {
        _output[0][_outOffset] = _input[0][_inOffset];
        _output[1][_outOffset] = _input[1][_inOffset];
      }
    } else {
      int halfSize = _size / 2;
      int nextStride = _stride * 2;
      ditfft2(_input, _inOffset, _output, _outOffset, halfSize, nextStride, _inversion);
      ditfft2(_input, _inOffset + _stride, _output, _outOffset + halfSize, halfSize, nextStride, _inversion);
      for (int i = 0; i < halfSize; i++) {
        int io1 = _outOffset + i;
        int io2 = io1 + halfSize;
        double real1 = _output[0][io1];
        double imag1 = _output[1][io1];
        double real2 = _output[0][io2];
        double imag2 = _output[1][io2];
        int iGonio = i * _stride;
        double sin = _inversion * sins[iGonio];
        double cos = coss[iGonio];
        double real3 = real2 * cos - imag2 * sin;
        double imag3 = real2 * sin + imag2 * cos;
        _output[0][io1] = real1 + real3;
        _output[1][io1] = imag1 + imag3;
        _output[0][io2] = real1 - real3;
        _output[1][io2] = imag1 - imag3;
      }
    }
  }

}
